#include "split_algorithm.hh"

static inline
mln::point2d
iterate_in_dir(mln::point2d& p,
	       e_dir dir)
{
  int dx = 0;
  int dy = 0;

  switch (dir)
    {
    case DIR_RIGHT:
      dx = 1;
      break;
    case DIR_RIGHT_UP:
      dx = 1;
      dy = -1;
      break;
    case DIR_UP:
      dy = -1;
      break;
    case DIR_LEFT_UP:
      dx = dy = -1;
      break;
    case DIR_LEFT:
      dx = -1;
      break;
    case DIR_LEFT_DOWN:
      dx = -1;
      dy = 1;
      break;
    case DIR_DOWN:
      dy = 1;
      break;
    case DIR_RIGHT_DOWN:
      dx = dy = 1;
      break;
    default:
      break;
    }

  return mln::point2d(p[0] + dx, p[1] + dy);
}

static inline e_dir inverse_dir(e_dir d)
{
  int a = (int)d;
  a = (a - 4 + 8) % 8;
  return (e_dir)a;
}

static inline
void
dec_dir(e_dir* d)
{
  int a = (int)*d;
  a = (a - 1) % 8;
  a = (a + 8) % 8;
  *d = (e_dir)a;
}

static inline
void
inc_dir(e_dir* d)
{
  int a = (int)*d;
  a = (a + 1) % 8;
  a = (a + 8) % 8;
  *d = (e_dir)a;
}

template <typename I>
static inline
e_dir
find_first_dir(I& ima,
	       mln::point2d& p)
{
  e_dir dir = DIR_NONE;

  mln::point2d n;
  do
    {
      dec_dir(&dir);
      n = iterate_in_dir(p, dir);
    }
  while (ima(n) != mln::literal::zero);

  while (ima(n) == mln::literal::zero)
    {
      inc_dir(&dir);
      n = iterate_in_dir(p, dir);
    }

  return (dir);
}

/*
  Next point in the contour, found either going clockwise or not,
  depending on the kind of iterator with which the function
  is templated with.
*/
template <typename I>
static
mln::point2d
next_pt_in_contour(I& ima,
		   mln::point2d& ori,
		   e_dir* last_dir)
{
  e_dir dir = (*last_dir);
  dec_dir(&dir);
  mln::point2d cur;

  do
    {
      cur = iterate_in_dir(ori, dir);
      dec_dir(&dir);
    }
  while (ima(cur) == mln::literal::zero);

  *last_dir = inverse_dir(dir);
  return cur;
}

/*
  The farthest point from line (begin, end) or NULL
  if that distance is inferior to the given precision.
*/
template <typename I>
static inline
mln::point2d
find_farthest(I& ima,
	      mln::point2d& begin,
	      mln::point2d& end,
	      float precision)
{
  float d = 0.;
  float dmax = 0.;
  e_dir dir = find_first_dir(ima, begin);

  float a = end[1] - begin[1];
  float b = begin[0] - end[0];
  float c = begin[1] * end[0] - begin[0] * end[1];
  float norm = sqrt ((a * a) + (b * b));

  mln::point2d max;
  mln::point2d cur;
  mln::point2d old;

  cur = begin;
  max = cur;

  while (cur != end)
    {
      old = cur;
      cur = next_pt_in_contour(ima, cur, &dir);
      
      d = a * cur[0] + cur[1] * b + c;

      if (d < 0)
	d = -d;

      if (d > dmax)
	{
	  dmax = d;
	  max = cur;
	}

    }

  if (dmax > precision * norm)
    return max;
  else
    return begin;
}

template <typename I>
void
split_rec(I& ima,
	  mln::point2d& begin,
	  mln::point2d& end,
	  std::list<mln::point2d>& l,
	  float precision)
{
  mln::point2d node = find_farthest(ima, begin, end, precision);

  if (node != begin)
    {
      split_rec(ima, begin, node, l, precision);

      l.push_back(node);

      split_rec(ima, node, end, l, precision);
    }
}

template <typename I>
void
split(I& ima,
      mln::point2d& begin,
      mln::point2d& end,
      std::list<mln::point2d>& l,
      float precision)
{
  mln::point2d node;

  node = find_farthest(ima, begin, end, precision);

  l.push_back(begin);

  if (node != begin)
    {
      split_rec(ima, begin, node, l, precision);

      l.push_back(node);

      split_rec(ima, node, end, l, precision);
    }

  l.push_back(end);

  node = find_farthest(ima, end, begin, precision);

  if (node != end)
    {
      split_rec(ima, end, node, l, precision);

      l.push_back(node);

      split_rec(ima, node, begin, l, precision);
    }
}

static inline
void
draw_polygon(mln::image2d<mln::value::rgb8>& ima,
	     std::list<mln::point2d>& l)
{
  if (l.size() == 0)
    return;

  mln::point2d prec = l.back();

  for (std::list<mln::point2d>::const_iterator i = l.begin();
       i != l.end(); ++i)
    {
      mln::draw::line(ima,
		      *i,
		      prec,
		      mln::literal::red);
      prec = *i;
    }
}

template <typename I>
static inline
std::pair<mln::point2d, mln::point2d>
get_initials(I& ima)
{
  mln::point2d begin(0, ima.ncols() / 2);
  mln::point2d end(ima.nrows() - 1, ima.ncols() / 2);

  while (ima(begin) == mln::literal::zero)
    begin[0]++;

  while (ima(end) == mln::literal::zero)
    end[0]--;

  return std::pair<mln::point2d, mln::point2d> (begin, end);
}

int
main (int argc, char** argv)
{
  mln::image2d<bool> ima;
  std::list<mln::point2d>* l;

  if (argc != 3)
    return 1;

  l = new std::list<mln::point2d> ();

  mln::io::pbm::load (ima, argv[1]);

  mln::image2d<bool>::piter p (ima.domain());

  std::pair<mln::point2d, mln::point2d> initials = get_initials(ima);

  split(ima, initials.first, initials.second, *l, atoi(argv[2]));

  mln::value::rgb8 v;
  mln::image2d<mln::value::rgb8> res = mln::data::convert (v, ima);

  draw_polygon(res, *l);

  mln::io::ppm::save(res, "out_"
		     + std::string(argv[2])
		     + std::string("_")
		     + std::string(argv[1]));

  return 0;
}
